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Bulk water presents a large number of crystalline and amorphous ices. Hydrophobic nanoconfinement is known to 
affect the tendency of water to form ice and to reduce the melting temperature. However, a systematic study of the ice 
phases in nanoconfinement is hampered by the computational cost of simulations at very low temperatures. Here we 
develop a coarse-grained model for a water monolayer in hydrophobic nanoconfinement and study the formation of ice by 
Mote Carlo simulations. We find two ice phases: low-density-crystal ice at low pressure and high-density hexatic ice at 
high pressure, an intermediate phase between liquid and high-density-crystal ice. 



1 Introduction 



1.1 The Model 
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Water phase diagram is very complex when compared 
to other liquids. For example, water is a polymorph with 
an unusually large number of solid phases (crystalline and 
amorphous ices): more than 20 with the last phase discov- 
erered in 2009 fT2\. Formation of ice in confined systems 
is a relevant subject in nanocience and biology, in areas like 
cryopreservation of food and human tissues or cells. Due to 
the property of water to expand when the liquid transforms 
into ice, the formation of ice in confinement can drammat- 
ically damage or destroy the confining structure. There- 
fore, it is important to understand the properties of ice in 
nanoconfinement, especially for hydrated systems at low 
temperature where water is highly confined, as for exam- 
ple in biological cells or on the surface of proteins at low 
hydration level ||71|9]. 

Simulations can help to answer open questions in this 
fields, but they are hampered by the large computational 
costs of calculations for large systems at low temperature 
with detailed models of water. With the aim of developing a 
model that preserve essential properties of water and is also 
computationally efficient, here we perform Monte Carlo 
simulations of a coarse-grained model of water. In its orig- 
inal formulation the model allows for very efficient simula- 
tion studies. On the other hand, it has a simple hamiltonian 
that allows for theoretical studies [5, 15 1, but does not al- 
low the study of structural properties because the positions 
of water molecules are coarse-grained. Here, we extend 
the model, introducing the coordinates of the molecules r,- 
to perform the structural analysis. We study the phase dia- 
gram and, in particular, how the hydrophobic nanoconfine- 
ment affects the ice formation for a water monolayer. We 
find two forms of ice and we characterize their structure. 

This work is organized as follows: we introduce the 
Model in the first section and give details about the Monte 
Carlo method in the second section; we present the results 
in the third section, discussing the structural and dynamical 
properties; we make our conclusions in the final section. 



We consider a water monolayer confined between two 
hydrophobic plates at separation h ~ Inm. The hydropho- 
bic interaction with the plates is schematically represented 
as purely repulsive. By molecular dynamics simulations of 
a detailed model of water, it has been shown that a water 
monloayer under these conditions forms a two-dimensional 
ice with square symmetry [8, 19, 18 1. Therefore, we adopt 
a square partition to coarse grain the confined water, di- 
vinding the total occupied volume V into cells of square 
section and height li. Each cell has a volume v = V /N and 
the coarse grain is made with the hypothesis that the sys- 
tem is homogenous and each cell contains one single water 
molecule. 

We consider the case in which pressure P, temperature 
T and number of molecules are fixed, and the total vol- 
ume V can vary. Therefore, the volume per molecule v, 
and the number density p = 1 / v, are functions of P and T 
at any fixed A^. 

The hamiltonian of the model has several water-water 
interaction terms. The first is the isotropic Van der Waals 
interaction, due to dispersive attractive forces and short- 
range repulsive interactions, represented by a Lennard- 
Jones potential: 
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(1) 



Here e = 5.8 kJ/mol is the characteristic energy, ro = 2.9 A 
is the diameter of the molecules and r,y the distance be- 
tween two water molecules in the cells / and j. In order 
to reduce the computational cost of the simulations, we in- 
troduce a cutoff distance at rcutoff = 2.5ro and add a linear 
term that set to zero the potential at rcutoff- 

Two neighbouring molecules can form a hydro- 
gen bond when the OH — O distance is less than 
rmax - roH = 3.14A, and if OOH < 30°. To account for 
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this interaction the model includes a term 

MiB - -JNhb (2) 

where J = 2.9 kJ/mol and A^hb = ^{(.;) A/ with 
Pij = Saij.aji&inj- rmax), and 0(x) = 1 if X > 0, or 
otherwise. Each molecules / has a bond indices 
cJij G {0,1,...^ — 1} for each nearest neighbor molecule 
j. The choice q = 180°/30° = 6 accounts correctly for 
the entropy loss associated with the formation of a hydro- 
gen bond because by definition Str^^.o^, = 1 if CJij — CJ,,, 
dojj.ajj = Otherwise. The notation (/, j) denotes that the 
sum is performed over nearest neighbors, implying that 
each molecule cannot form more than 4 bonds with its near- 
est neighbours. 

When water molecules form a hydrogen bond network, 
the resulting configuration occupies more space than at 
close packing. This effect is included in the model as a 
volume increase per formed bond equal to vhb ~ O.Svq, 
correspondig to the average density increase between high 
density ices VI and VIII and low density (tetrahedral) ice 
Ih in bulk water 1(131 1741 . The total volume of the system is, 
therefore, V —Vo + vhbMib, where Vq is the volume when 
there are no hydrogen bonds. 

As an effect of cooperativity, the O-O-O angle distri- 
bution becomes sharper at lower T, reducing the possible 
orientations of the molecules IfTTl |3] |6] . This cooperative 
term, resulting from three-body interactions, is accounted 
for by the term 

'■ ik-l), 

where Ja = 0.29 kJ/mol, and {k, £)i indicates each of the six 
different pairs of the four bond-indices a,; of a molecule 
/. The effect of this term is to locally drive the molecules 
toward an ordered configuration. 

In its original formulation the model is defined by 
coarse-graining the molecules coordinates r,y with the cen- 
ter of each cell. Furthermore, the effect of the cooperativ- 
ity on the O-O-O angle distribution is taken into account in 
terms of the associated entrophy change, but not in terms 
of angular coordinates. Therefore, no detailed structural 
analysis is possible. To allow the calculation of the radial 
distribution function g{r) and the angular distribution func- 
tion g{9), in the following subsection we extend the orig- 
inal model introducing a term that explicitly depends on 
these variables. 

1.1.1 Extensionofthemodel 

The new Hamiltonian term is a three-body interaction 
that depends on the formation of hydrogen bonds between 
triads of molecules and their relative angles 0^,,: 

^Mk^Jel, L PiMei,), (4) 

where Jg = 0.5e. The sum is over all the neighbouring 
pairs of molecules k and / that are bonded to the molecule 



/, with the restriction that k and / must be second nearest 
neighbors to each other The function A(0) is a smooth 
function of the angle between the centers of the three 
molecules with a minimum at n/2. We adopt this choice 
because molecular dynamics simualtions of a detailed wa- 
ter model show that, under the conditions considered here, 
a confined water monolayer forms a square crystal ||8] [19] . 
We chose 

A(0) = i[l-cos(40)], 

which is a non-negative function in [0, 1] with minima at 
k/2 and n, and we approximate it with 

A(0) ~4(0-7r/2)2 (5) 

around ~ 7r/2. 

This value of Jg is set to avoid the formation of bonds 
when the « 60°, because 

.JfeiO'j, ^7t/3) = J% {9}, = 2;r/3) ~ 4NJe , 

and 

J^+J^j- ANJg + 2NJ = 2Ne-Ne=Ne>0. 

Therefore, the formation of hydrogen bonds is energeti- 
cally unfavourable when « 60°. 

The total hamiltonian of the model is 

Jif = + '^^B + + --^e ■ (6) 
1.2 Metropolis MC method 

We perform MC simulations at constant number of 
molecules = 900 and fixed pressure P and temperature 
T, allowing fluctuations of the volume V. One MC step 
consists in updating 5A^+ 1 variables: vectors r,- describ- 
ing the position of molecules with respect to the center of 
their cell, 4A^ bondig indices Gij and the total volume Vq. 
We adopt the Metropolis algorithm: we choose one of the 
5A^ + 1 variables at random and attempt to change its state 
to a new random value. We accept the new state with prob- 
ability exp[— )3AG] if AG > and with probability 1 oth- 
erwise. Here /3 = l/ksT, kg is the Boltzamann constant, 
AG = Gnew — Gold is the change in Gibbs free energy if the 
new state is accepted, and 

G = U +PV -TS 

^,yif + PV -NkBT\og{V). (7) 

For a bonding index Oij, the new state is choosen at 
random among the q — 6 possible states. For each com- 
ponents of rjj, the new value is set to ra.new = '"a.oid + £r, 
where Sr G [—dr,+5r] is a random number and a — x, y 
(we do not change the component z and consider it as a 
coarse-grained variable). The volume is updated with a 
random change V,,"'^* Vg"''* + Sy where Sy G [-5V, +5V] 
is a random number [10|. The parameters 5r and dV are 
adapted in such a way to keep the acceptance ratio « 40% 
(adaptive step size algorithm) Iil6t irj. 
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At any P, we equilibrate the system from random con- 
figurations at high T for 10^ ^ 10^ MC steps and calculate 
the thermodynamic averages over the following 10*^-10^ 
MC steps. Keeping P constant, we perform an annealing, 
i.e. we decrease the temperature T a few K and, starting 
from the last configuration at the previous temperature, we 
use the same statistics for equilibration and calculation of 
the averages. To take into account the correlation of the 
data for the calculation of the error on the estimates, we 
perform blocking averages where the size of each block 
depends on P and T and is determined as twice the num- 
ber T of MC steps needed to have uncorrected data. The 
number T is estimated from the autocorrelation functions 
introduced in Section|4] 



2 Results 

2.1 Phase Diagram 

The phase diagram of the model displays the liquid-gas 
first-order phase transition ending in a critical point (FiglTJ. 
In the liquid phase we observe that at any P < 0.2 GPa the 
density is non monotonic. The locus of temperatures of 
maximum density (TMD) follows a line in the P-T phase 
diagram, reproducing one of the characteristc anomalies of 
water 

At low P and low T we find a rapid decrease of den- 
sity p. This is the consequence of formation of a large 
number of hydrogen bonds and the cooperative reorienta- 
tion of the molecules into a crystal configuration. In the 
next section we caracterize this crystal as low-density crys- 
tal (LDC) with a square cell in its 2D projection. 

At high P and low T our structural and dynamical anal- 
ysis, presented in the next sections, show that the system 
"freezes" into a solid. However, the solid has no long- 
range translational order, but short-range translational or- 
der and quasi-long-range orientational order This is, by 
definition, an "hexatic" phase, described by the theory of 
Kosterlitz, Thouless, Halperin, Nelson, and Young |2 11 1 
for crystallization in 2D systems. The theory tells us that 
the hexatic phase is intermediate between the crystal and 
the liquid phases and is separated by continuous phases 
transitions with both phases. This is consistent with the 
fact that we do not observe any discontinuity in the den- 
sity at high P and low T (Fig. [U and that our system is 
essentially in 2D because we coarse-grain the z-component 
of the molecules. The hexatic-liquid coexistence is charac- 
terized by the unbinding of disclinations, i.e. lines of de- 
fects at which rotational symmetry is violated. The crystal- 
hexatic coexistence is where occurs the unbinding of dis- 
locations, i.e. particle-like topological defects. The asso- 
ciated crystal phase is characterized by the same orienta- 
tional order of the hexatic phase, that is, as described in the 
next section, hexagonal (or close-packing) and has a higher 
density of the LDC. We therefore call it high-density crys- 
tal (HDC). Finally, the structural analysis allows us to esti- 
mate the coexistence line between the LDC and the hexatic 
phase and the triple point where liquid, hexatic and LDC 
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Figure 1 : (a) Isobars for the confined water monolayer in 
the p — r phase diagram display a discontinuity that marks 
the first-order liquid-gas phase transition ending in a crit- 
ical point (at P = 0.16 ±0.02 GPa, T = 1500 ±200 K) 
and a density maximum for P < 0.20 GPa. Isobars are for 
P e [0.02,0.3] GPa separated by 0.04 GPa (from the bot- 
tom), (b) At low T a sudden decrease in p suggests the 
occurrence of another phase transition. As we discuss in 
the analysis, we associate this sudden change in p to the 
first-order phase transition between the liquid and a low 
density crystal (LDC). Pressures are as in (a) but with a 
separation of 0.02 GPa. (c) The P —T phase diagram with 
gas-liquid coexistence line ending at the gas-liquid critical 
point (open circle) and the line of temperatures of maxi- 
mum density (TMD). At low T and low P, the liquid-LDC 
coexistence line merges, at about f = 0. 11 ± 0.01 GPa and 
r = 190±20 K (open triangle), with a line with positive 
slope corresponding to the liquid-hexatic phase transition, 
as described in the text. Full circles are the state points of 
the distribution functions in Fig.s l2l3l 
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phases coexist (Fig. [T]i. 

It is interesting to observe that the phase diagram of 
our confined monolayer reproduces the change of slope of 
the "ice" line observed for bulk water. The slope is neg- 
ative at low P and is positive at high P. The ice phase at 
low P is LDC characterized by hydrogen bonds at 90°. At 
high P the solid phase is hexatic, where the number of hy- 
drogen bonds is largely reduced and the water interaction 
is dominated by the Lennard-Jones potential, as in simple 
liquids. 

2.2 Structural Properties 

The radial distribution function 

To study the static properties of the system we calcu- 
late the radial distribution function (RDF) as 



(8) 



The quantity g{r)2nrdr is proportional to the probability 
of finding a molecule at a distance r from a central one. 

By crossing the ice-liquid lines of FigH]; we can ob- 
serve structural changes in the g{r). At low P (Fig|2^), 
we find a large change in g{r) within a narrow range of T, 
marking the occurrence of the liquid-LDC first-order phase 
transition. The LDC is characyerized by square long-range 
translational order. 

At high P (Fig|2]3), we observe a shoulder in the sec- 
ond peak of the g{r) for the liquid. This shoulder develops 
into a small peak at lower T . This structural change has 
been characterized fTT] as the liquid-hexatic second-order 
phase transition. This interpretation is consistent with the 
analysis of the typical configurations at the lower T, show- 
ing liquid-like short-range translational order and crystal- 
like long-range orientational (hexagonal) order The hex- 
atic phase is the precursor of the HDC close-packing crys- 
tal. The solid-like properties of this phase are confirmed by 
the analysis presented in the next section. 

Finally, at low T by increasing P the structural analys 
allows us to locate the coexistence between the LDC and 
the hexatic phase. The transition is charcterized by a sharp 
change of g{r) indicating a first-order phase transition be- 
tween the LDC and the hexatic phase. 

The angular distribution 

To charcterize the structure, we also calculate the 
O-O-O (not normahzed) angular distribution function, de- 
fined as 



O0('-«-'-max)5(0-0 



g{B)^'-Y^Y^&{r,- 

(9) 

The quantity g{9)r dr dO is proportional to the probability 
of finding two molecules (j, k) at a distance r < rmax from 
a central one / and forming an angle = 0. The condi- 
tion r < Tmax limit our calculation to the first shell, in the 
condensed phase, of the central molecule. 
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Figure 2: Radial Distribution Functions for the state 
points marked as coloured full circles in Fig{T] (a) At 
P = 0.06 GPa and T = 246 K (in the Hquid phase) and 
r = 218 K (in the LDC phase). Here and in the next panel, 
insets show a portion of typical configurations at the state 
points represented in the main panels. The g{r) for the 
LDC has many peaks correspoding to the long-range trans- 
lational order of the square crystal, while the g{r) for the 
liquid near the coexistence shows precursors of the LDC 
stmcture. (b) At P = 0.24 GPa and T = 274 K (in the liq- 
uid phase) and T = 218 K (in the hexatic phase). The liq- 
uid g{r) has a shoulder in the second peak that splits into 
a small peak in the hexatic phase. The hexatic phase has 
liquid-hke short-range translational order due to the pres- 
ence of many disclinations, but crystal-like long-range ori- 
entational order, emphasized by links in the inset describ- 
ing the hexatic phase. 
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Figure 3: The angular distribution functions g{9), calcu- 
lated at the same state points as in Fig|2] emphasizes the 
long-range orientational order in the LDC and the hexatic 
phases. 

Our analysis of g{9) (Figl3]l is complementary to that 
of g{r) and emphasizes the appearence of the long-range 
orientational order in the LDC and the hexatic phases. In 
the two solid phases, the positions of the peaks are related 
to the symmetry of each crystal phase: the square LDC 
structure has peaks at 90° and 180°, and the peaks of the 
hexatic solid phase, with the same symmetry as the HDC, 
are centered at 60°, 120° and 180°. In the liquid phase, 
instead, g{9) never goes to zero showing the absence of 
orientational order 

2.3 Dynamical Properties 

The study of the autocorrelation functions provide rel- 
evant informations about both the MC dynamics and the 
transport properties of the system. From them we extract 
the correlation times necessary to calculate in a correct way 
the statistical errors of our observables. Moreover, from the 
correlation time we estimate when the dynamics of the sys- 
tem is liquid-like or solid-like. 

We calculate the hydrogen bonds autocorrelation func- 
tion 



Cm(0 



]_y{Mi{tO 



-t)Mi{to))-{M,Y 



N 



(Mf) - {Mi] 



(10) 



of the average molecular bonding index of molecule / 
Mi = J Y.j Oij — {q—l)/2. This quantity describes the hy- 
drogen bonds dynamics of water molecules. 

Next, we calculate the translational autocorrelation 
function 



Cr{t) 



}_y^h{t0 + t)-5ri{tQ))-{5r,f 

AT i—i 



N 



{5r})-{5hY 



(11) 



where 5r = r, — tq , is the displacement of each molecule / 
from the center of its cell. In Eq. (fTOt and (fTTT i the time t is 
measured in MC steps and can be related to real time only 
by comparison with experiments. For example, it can be 
shown that the conversion factor between a MC step and 
real time unit rescales logarithmically with T at ambient 
pressure |i9J. 
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(a) HB autocoiTelation 
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Figure 4: Autocorrelation functions of the average hy- 
drogen bond index M (a), and of the displacement of a 
molecule with respect to the center of its cell (b). The func- 
tions are calculated at the state points indicated in the leg- 
end (with T expressed in K and P in GPa). Pressures are 
as in Fig. 12] and, for both pressures, the water monolayer 
is Hquid for T > 218 K and solid for T < 218 K: LDC at 
P = 0.06 GPa and hexatic at P = 0.24 GPa. 



For each quantity we define the correlation time T 
as the time at which the normalized correlation function 
Eq. ( [Tol l and (fTTI) decay to 1 /e. Our calculations show 
that at low P, the hydrogen bond correlation function Cm 
(FigH^) is exponential in the liquid phase, but has a non- 
exponential behavior in the LDC phase, with a correlation 
time T that largely increases for decreasing T, as expected 
in the solid phase characterized by a well developed hy- 
drogen bond network [4]. By increaing P, the number of 
hydrogen bonds largely decreases and the hydrogen bond 
correlation function Cm shows a much faster decay to zero. 
Nevertheless, at high P and low T the function is not- 
exponential consistent with the approach of a frozen state. 

For the translational autocorrelation function Cr 
(Fig|4j3) we find a much slower decay to zero at all the 
state points. For the state points corresponding to the solid 
phases, at low T and any P, the function has an evident 
non-exponential behavior and an extremely long correla- 
tion time T, two orders of magnitude greater than Cm- This 
is consistent with the arrested translational dynamics of the 
solid phases. 

Next, we analyze the behavior of the variance (5^M,) 
of Mi, and (5^(5r,)) of 5r,, defined as the normalization 
factors of Eq.s ( fTOl i and (fTTl i. respectively (Fig|5]i. 

As expected, at low P both variances have discontinu- 
ties at the temperature of the liquid-LDC first-order phase 
transition. The increase of (5^M,) for decreasing T is the 
consequence of the increase of the fluctuations in the hy- 
drogen bonds network at low T . 

On the other hand, the translational variance decreases 
for decreasing T and presents discontinuities at the first- 
order phase transitons, e. g. being crystal at T < 274 K, 
gas at r > 483 K and a liquid at the intermediate values 
at P = 0.02 GPa. At high P, instead, the slowing down of 
the translational dynamics occurring at the liquid-hexatic 
coexistence is marked by a non monotonic (5^(5r,)), sug- 
gesting an out of of equilibrium behavior at very low T . 
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Figure 5: The variance of the average hydrogen bond index 
Mi at (a), and of the displacement of the molecule in- 
side the cell at (b). At low P both display a discontinuous 
decrease at the liquid-LDC first-order phase transition. At 
high P, around the liquid-hexatic second-order phase tran- 
sition, (S^Mj) increases continuously, while (5^(5r,)) has 
a discontinuity. At high T, (5^(5r,)) shows the discontinu- 
ity associated to the liquid-gas first-order phase transition. 
Pressures are expressed in GPa. 

3 Conclusions 

We study by efficient Monte Carlo simulations a 
coarse-grained model for a water monolayer in hydropho- 
bic nanoconfinement and find two forms of ice at low T. 
At low pressure, the model reproduces the occurrence of 
low-density-crystal (LDC) ice with hydrogen bonds form- 
ing a square network, as observed with detailed molecu- 
lar dynamics simulations. At high pressure, where detailed 
molecular dynamics simulations are not available, we find 
a hexatic ice, separated from the liquid phase by a second- 
order phase transition. Our structural analysis shows that 
the hexatic phase has solid-like long-range orientational or- 
der and liquid-like short-range translational order 

By studing the autocorrelation functions and the vari- 
ances of bonding and translational parameters, we observe 
different behaviors at low and high P that we can relate 
to the thermodynamic phases. The dynamics at low T be- 
come slow at any pressure, and possibly fall out of equi- 
librium at low T in the hexatic phase, before entering the 
high-density-crystal (HDC) ice phase. 
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